De novo emergence of a remdesivir resistance mutation during treatment of persistent SARS-CoV-2 infection in an immunocompromised patient: a case report

SARS-CoV-2 remdesivir resistance mutations have been generated in vitro but have not been reported in patients receiving treatment with the antiviral agent. We present a case of an immunocompromised patient with acquired B-cell deficiency who developed an indolent, protracted course of SARS-CoV-2 infection. Remdesivir therapy alleviated symptoms and produced a transient virologic response, but her course was complicated by recrudescence of high-grade viral shedding. Whole genome sequencing identified a mutation, E802D, in the nsp12 RNA-dependent RNA polymerase, which was not present in pre-treatment specimens. In vitro experiments demonstrated that the mutation conferred a ~6-fold increase in remdesivir IC50 but resulted in a fitness cost in the absence of remdesivir. Sustained clinical and virologic response was achieved after treatment with casirivimab-imdevimab. Although the fitness cost observed in vitro may limit the risk posed by E802D, this case illustrates the importance of monitoring for remdesivir resistance and the potential benefit of combinatorial therapies in immunocompromised patients with SARS-CoV-2 infection.

Clinical and virological response to remdesivir and mAb therapy. The patient was hospitalized at our institution and received ten-day course of RDV (day 148-157), which was extended from the standard five-day regimen, because of concerns for clinical and virologic relapse in immunocompromised hosts with persistent infection. The patient had a clinical response to therapy, characterized by resolution of her fever (day 149, Fig. 1a), normalization of CRP (day 156, Fig. 1c) and improvement in the opacities on chest CT (day 162, Fig. 1d). She also had an initial virologic response to therapy with an increase in RT-PCR Ct values, decreasing numbers of plaque forming units (PFU) in viral cultures and lower proportions of SARS-CoV-2 subgenomic RNA (sgRNA), as measured in amplicon sequencing data, in her respiratory tract specimens by day 152 (Fig. 2a).
However, a recrudescence of viral shedding occurred during and after RDV therapy, as evidenced by decreasing RT-PCR Ct values which reached 18 on day 160 and increasing PFU and sgRNA in respiratory tract specimens on day 156 (Fig. 2a). The patient received an 8 g infusion of casirivimab-imdevimab 6,7 on day 163 after expanded access for the use of the monoclonal antibody (mAb) therapy was approved. A rapid and sustained virologic response was observed with undetectable cultured virus, low sgRNA and negative NP RNA in specimens by days 164, 166 and 217, respectively. The patient's anosmia resolved roughly 17 days after administration of casirivimab-imdevimab (day 180. Fig. 1a). During the following five-month period (days 160-292), the patient did not have recurrence of COVID-19 related symptoms and had a chest CT which showed minimal opacities (day 170, Fig. 1d), and had high anti-SARS-CoV-2 S1 protein IgG titers (Fig. 2a). Peripheral blood neutrophil counts, hemoglobin levels and serum inflammatory markers normalized during convalescence (Fig. 1c).
Emergence of the nsp12 E802D mutation during remdesivir therapy. To identify SARS-CoV-2 mutations that arose during the course of illness, we longitudinally sampled multiple sources of the patient's tissues and secretions (Supplementary Table 1) and performed WGS on an Illumina NextSeq platform using the Swift SARS-CoV-2 multiplex amplicon sequencing panel 8 . Viral N1 or N2 RNA was detected in 19 of 21 nasopharyngeal, 10 of 10 saliva, 10 of 10 stool and 0 of 10 whole blood specimens. Among the 39 samples with detectable N1 or N2 RNA, we assembled 27 whole SARS-CoV-2 genomes from nasopharyngeal (12), saliva (9) and stool (6) specimens. One of the nasopharyngeal specimens was obtained from the initial phase of illness (day 36, Fig. 2a). Phylogenetic analysis found that sequenced genomes belonged to a single lineage (Pango B1) within Nextstrain clade 20 C ( Supplementary Fig. 1), indicating that viral genomes identified during the course of illness were derived from intra-host diversification following infection with a single strain.
Analysis of assembled viral genomes identified a mutation, E802D, in nsp12, whose detection in patient specimens was temporally associated with RDV therapy (Fig. 2b). E802D was not identified at an allele frequency above 1% in either the specimen obtained during the initial phase of illness (day 36) or in specimens collected during the first 5 days (days 148-152) of RDV therapy. The mutation was first detected 7 days after initiation of RDV therapy (day 155) and accounted for 23% and 96% of the allele frequency in nasopharyngeal and saliva specimens, respectively, by day 160. In a PAML analysis 9 of the Orf1a/1ab reading frame, this position was identified as having the highest probability (BEB Prob(w > 1), 0.984; mean omega, 9.90) for positive selection after RDV administration.
The E802 residue of nsp12 resides in the palm sub-domain which contains a portion of the residues that comprise the active site of the SARS-CoV-2 RNA polymerase 1,2,10 . The E802 residue participates in an electrostatic network with D804 and K807 which stabilizes the loop involved with binding to the nascent RNA (Fig. 2c, Supplementary Fig. 2a). E802D has been identified in an in vitro RDV resistance selection experiment and was found to confer a~2.5-fold increase in IC 50 to the drug 2 . Together, the temporality of E802D emergence in the patient, its location in nsp12 and the in vitro identification of the same mutation with a RDV resistance phenotype support the plausibility that that RDV treatment of the patient selected for variants with the E802D mutation, which in turn contributed to the observed rebound viral shedding.
In vitro validation of nsp12 E802D as a remdesivir resistance mutation. To validate the resistance phenotype of E802D, we engineered this mutation and an E802A mutation, which had been also shown to confer RDV resistance 2 , into an infectious molecular clone of SARS-CoV-2/WA01 (icSARS-CoV-2-mNG), which expresses the mNeon Green reporter and is ORF7a depleted 11 . Replication kinetics of the E802D and E802A mutants in Vero-E6 cells revealed decreased viral replication relative to parental icSARS-CoV-2-mNG (Fig. 2d), suggesting that these mutations may impart a fitness cost. In the presence of high concentrations of RDV (5 µM), the E802D mutant replicated to higher titers than the parental virus ( Supplementary Fig. 2b), confirming that this mutation confers resistance to RDV. RDV dose-response curves demonstrated that E802D and E802A mutants had significant increases in IC 50 values (4.2 μM and 2.7 μM, respectively), relative to parental icSARS-CoV-2-mNG (0.7 μM) (Fig. 2e). Assessment of RDV cytotoxicity demonstrated limited impact on cell death rates below 10 μM, indicating that the RDV resistance phenotype conferred by mutations at E802 may occur at physiologically relevant concentrations (Supplementary Fig. 2c).
These findings suggest that RDV resistance mediated by substitutions at residue 802 in nsp12 may lead to therapeutic failure in the setting of an immunocompromisted host. We speculate that mutations at 802 distort the active site in a way that either enables the enzyme to exclude RDV or alleviates the steric clash mediated by S861, thereby allowing the enzyme to escape RDV-mediated chain termination 1,10 . Further biochemical and structural characterization will be needed to delineate the mechanism by which substitutions at residue E802 confer resistance to RDV.
The therapeutic impact of RDV resistance mutations at E802 may be attenuated by a fitness cost imparted by these mutations, which we observed in vitro. Although RDV therapy has been widely administered to patients during the pandemic, E802D and all substitutions at residue 802 have been found in 131 and 297, respectively, of the 4.8 M genome sequences obtained from patient isolates in the GISAID database (www.gisaid.org; accessed 11/07/21). Consistent with the observation of a fitness cost, the allele frequency of E802D in nasopharyngeal specimens from the patient decreased from 22.6% to <1% between the period (day 160-162) after completion of RDV therapy and prior to illness which include the period when remdesivir (gray shading) and casirivimab/imdevimab (light blue shading) were administered. b Allele frequencies of E802D in nsp12 as ascertained by whole genome sequencing. c Crystal structure (PDB:7BV2) depicting nsp12 E802 (dark blue), its putative hydrogen bonds (yellow dashed lines) within the palm domain (yellow), and the residues (green) that interact with the replicating RNA molecule (orange). Remdesivir monophosphate (RMP) is depicted in light blue. d Viral growth curves from icSARS-CoV-2 mNG WT, E802D (patient), and E802A (control) nsp12 mutants on a ORF7a depleted backbone. Results are depicted as mean and standard deviation (error bars) of data from biological replicates (n = 3). Significance between WT and E802D mutants was assessed by unpaired, two-sample t tests, **p < 0.01 (p = 0.002 at 24 h, p = 0.006 at 48 h). e SARS-CoV-2 inhibition by RDV assessed at 48 h post-infection (0.01 multiplicity of infection) as determined by quantitative image analysis of percentage of cells expressing mNG. Color-coded curves represent a non-linear least squares fit of biological replicates (n = 3) and shading represents 95% confidence intervals of the fit. casirivimab-imdevimab administration (Fig. 2b). However, E802 allele frequencies were significantly higher (50-96% at days 158-167) in saliva specimens during this period and did not decline until casirivimab-imdevimab was administered.
An explanation for these discordant findings is unclear and may relate to the differential compartmentalization of variants within the respiratory tract or the stochastic nature of sampling the nasopharyngeal compartment as opposed to saliva where virus from a wider geographic distribution may pool. A sufficient number of nasopharyngeal and saliva specimens were not obtained to perform a compartmentalization analysis. Given that saliva is infrequently sampled, RDV resistance mediated by E802 mutations may be under-detected. Furthermore, analysis of the 131 genomic sequences that have been found in GISAID to contain the E802D mutation identified two instances where sequences were geographically and temporally clustered, suggesting that local transmission of E802D variants may have occurred ( Supplementary Fig. 3). Continued surveillance will thus be needed to elaborate whether the occurrence of E802D, as well as other substitutions at this residue, will be limited because of an underlying fitness cost or alternatively, emerge to pose a broader risk for RDV resistance.
Additional mutations that were temporally associated with remdesivir therapy. We did not identify nsp12 mutations at sites other than residue 802 which were temporally associated with RDV therapy. A mutation at D484Y has been reported in a posttreatment sample from an immunocompromised patient who failed RDV therapy 12 , but the RDV resistance phenotype associated with this mutation, whose location is distal to the active site, has not been confirmed. We identified two additional mutations, A504V in nsp14 (exonuclease) and I115L in nsp15 (endoRNAse), whose allele frequencies increased during and after RDV treatment ( Supplementary Fig. 4). However, A504V in nsp14 was present (allele frequency 85.1%) in a specimen obtained during her initial phase of illness (day 36). I115L in nsp15 was not detected in the early-phase specimen; we therefore cannot exclude the possibility that this mutation in the endoR-NAse gene may contribute to RDV resistance. In addition, further investigation is needed to evaluate whether epistatic interactions influence the expression of a remdesivir resistance phenotype beyond the isolated contribution of substitutions at E802D in nsp12.

Discussion
The case highlights several issues for the management of SARS-CoV-2 infection in immunocompromised patients. First, the course of infection was notable for being indolent despite a six-month period of sustained high-grade viral shedding and underlying lung involvement. We identified six mutations, including a frameshift mutation in ORF3a (N257fs), which were present in the consensus genomes throughout the course of illness prior to RDV treatment ( Supplementary Fig. 5). Yet none provides a clear explanation for an attenuated phenotype. Longitudinal profiling of PBMC demonstrated profound depletion of circulating B-cell populations and a T-cell phenotype consistent with chronic antigen exposure as evidenced by the increased activation (CD38+, HLA-DR+) and exhaustion (PD-1, Tim-3) markers ( Supplementary Fig. 6, gating strategy in Supplementary Fig. 8). Given reports that maladaptive immune responses contribute to the risk of severe COVID-19 13,14 , we speculate that the compromised immune response was unable to control viral load yet conversely, protected the patient from developing severe complications.
Second, SARS-CoV-2 infection was associated with the development of severe refractory neutropenia and moderate anemia in the patient (Fig. 1c). While we are unable to discern if the cytopenia was due to direct or indirect effects of infection, SARS-CoV-2 has been detected in the bone marrow of an immunocompromised patient who developed pancytopenia during infection 15 . Our patient's cytopenia resolved after casirivimabimdevimab therapy, suggesting a causal link between SARS-CoV-2 infection and cytopenia in this setting.
Third, the case adds to the evidence from prior studies 4,5,12 that monotherapy with RDV may be associated with a risk for clinical or virologic relapse in immunocompromised patients. Of note, the patient experienced a marked clinical response to RDV therapy, and we were unable to assess whether residual complications (i.e., anosmia, anemia, neutropenia) would have resolved without the intervention of mAb therapy. Nevertheless, the patient developed virologic failure during and after RDV therapy, which prevented her from terminating her prolonged isolation.
Finally, the case illustrates the potential benefit of mAb therapy for immunocompromised patients with persistent SARS-CoV-2 infection. Prior to initiation of mAb therapy, two nonsynonymous mutations, H655Y (Supplementary Fig. 7a) and P812L ( Supplementary Fig. 7b), emerged which were not present in the initial specimen (day 36), indicating that spike protein continues to evolve and adapt to the host in the absence of selective pressure from an adaptive immune response. One week after casirivimab-imdevimab was administered (day 170), we identified a new spike protein mutation, A348S in the RBD domain, which persisted to day 177 ( Supplementary Fig. 7c). WGS was not performed on specimens collected after day 177 because of low amounts or undetectable viral RNA. A348S is not proximal to the mAb binding sites (Supplementary Fig. 7d) 16 and was not identified in a comprehensive analysis of resistance mutations 6 . In the setting of a negligible cellular immune response, mAb therapy cleared viral shedding, abrogated the residual complaint of anosmia, and resolved the blood dyscrasias.
In summary, we identified the de novo emergence of a RDVresistance mutation, E802D, following initiation of RDV in an immunocompromised patient with persistent SARS-CoV-2 infection. While the finding is limited to a single case and requires confirmation of its generalizability in larger patient populations, it suggests that RDV may impart a selective pressure in vivo to drive evolution of the virus. E802D is associated with a fitness cost in vitro, which may limit the broader impact of this mutation on the development of secondary resistance during treatment and the risk for primary resistance through transmission of resistant variants. Yet, our findings underscore the importance of immunocompromised hosts with uncontrolled viral replication as a source of genetic diversification 3,4,17-19 and selection of mutations that may potentially impart adverse consequences for antiviral therapy. Enhanced genomic surveillance of immunocompromised patients may thus be warranted. As observed in this case, initiation of anti-SARS-CoV-2 mAb may serve as a therapeutic option to achieve rapid and sustained virologic responses and improved clinical outcomes in immunocompromised patients.

Methods
Ethics declaration. Regulatory approval to administer an 8 g dose of casirivimab/ imdevimab was obtained from the U.S. Food and Drug Administration under investigational new drug application no. 153390. This study was also approved by Yale Human Research Protection Program Institutional Review Boards (FWA00002571, protocol ID 2000027690, 2000029277). Informed consent was obtained from both the patient and healthcare workers in accordance with CARE guidelines and in compliance with the Declaration of Helsinki principles.
Isolation of patient plasma and PBMCs. Patient whole blood was collected in sodium heparin-coated vacutainers and gently agitating at room temperature until sample pick-up by IMPACT team members. All blood was processed on the day of NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29104-y ARTICLE NATURE COMMUNICATIONS | (2022) 13:1547 | https://doi.org/10.1038/s41467-022-29104-y | www.nature.com/naturecommunications collection. Plasma samples were collected after centrifugation of whole blood at 400 g for 10 min at room temperature (RT) without brake. The undiluted serum was then transferred to 15-ml polypropylene conical tubes, and aliquoted and stored at −80°C for subsequent analysis. PBMCs were isolated using Histopaque (Sigma-Aldrich, #10771-500 ML) density gradient centrifugation in a biosafety level 2+ facility. After isolation of undiluted serum, blood was diluted 1:1 in room temperature PBS, layered over Histopaque in a SepMate tube (StemCell Technologies; #85460) and centrifuged for 10 min at 1200 g. The PBMC layer was isolated according to the manufacturer's instructions. Cells were washed twice with PBS before counting. Pelleted cells were briefly treated with an ACK lysis buffer for 2 min and then counted. Percentage viability was estimated using standard Trypan blue staining and an automated cell counter (Thermo-Fisher, #AMQAX1000).
SARS-CoV-2 specific-antibody measurements. Triton X-100 and RNase A were added to serum samples at final concentrations of 0.5% and 0.5 mg/ml, respectively, and incubated at room temperature (RT) for 30 min as described 20 . 96-well MaxiSorp plates (Thermo Scientific #442404) were coated with 50 μl/well of recombinant SARS Cov-2 S1 protein (ACROBiosystems #S1N-C52H3-100ug) at a concentration of 2 μg/ml in PBS and were incubated overnight at 4°C. The coating buffer was removed, and plates were incubated for 1 h at RT with 200 μl of blocking solution (PBS with 0.1% Tween-20, 3% milk powder). Serum was diluted 1:50 in dilution solution (PBS with 0.1% Tween-20, 1% milk powder) and 100 μl of diluted serum was added for two hours at RT. Plates were washed three times with PBS-T (PBS with 0.1% Tween-20) and 50 μl of HRP anti-Human IgG Antibody (Gen-Script #A00166, 1:5,000) diluted in dilution solution added to each well. After 1 h of incubation at RT, plates were washed three times with PBS-T. Plates were developed with 100 μl of TMB Substrate Reagent Set (BD Biosciences #555214) and the reaction was stopped after 15 min by the addition of 2 N sulfuric acid. Plates were read at a wavelength of 450 nm and 570 nm.
Cell lines and virus. Vero E6 kidney epithelial cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 1% sodium pyruvate (NEAA) and 5% fetal bovine serum (FBS) at 37°C and 5% CO2. The cell line was obtained from the ATCC and has been tested negative for contamination with mycoplasma. SARS-CoV-2, strain USA-WA1/2020, was obtained from BEI Resources (#NR-52281) and was amplified in Vero E6 cells. As described 21 , cells were infected at a MOI 0.01 for four three days to generate a working stock and after incubation the supernatant was clarified by centrifugation (450 g × 5 min) and filtered through a 0.45-micron filter. The pelleted virus was then resuspended in PBS then aliquoted for storage at −80°C. Viral titers were measured by standard plaque assay using Vero E6 cells. All experiments were performed in a biosafety level 3 with the Yale Environmental Health and Safety office approval.
Virus titration from nasopharyngeal swabs. SARS-CoV-2 titers from longitudinal nasopharyngeal swabs were assessed by conventional plaque assay. Nasopharyngeal swabs were collected and swabs stored and frozen in universal viral transport media (VTM) prior to plaquing. Each nasopharyngeal sample was ten-fold serially diluted (dilution range from 1:10 to 1:1000000). Diluted NP samples were added on a monolayer of Vero E6 cells. After 1 h of infection, DMEM (2% FBS), and 1.2% of Avicel (Dupont, 9004346) were overlaid. 72 h post inoculation, media was removed and 4% Paraformaldehyde added for 1 h. Cells were stained with 0.1% of crystal violet to visualize SARS-CoV-2 plaques. Plaques were counted by single operator. Each sample was tested in duplicate.
RNA extraction and PCR. RNA was extracted on a Qiagen Bio Robot (Qiagen, Hilden, Germany) using the QIAamp Virus Bio Robot MDx Kit (Qiagen) following manufacturer's instructions. RT-PCR was performed using AgPath-ID One-Step RT-PCR kit (Life Technologies, Carlsbad, CA, USA) with CDC primers and probes for N1/N2 following protocol and concentrations defined previously 22 . Each 25 µL of reaction mix included 2X RT-PCR buffer, 25X enzyme mix, primers-probes, and 5 µL of the extracted RNA. Cycle parameters on the ABI 7500 real-time PCR system (Applied Biosystems) were 10 min at 48°C for reverse transcription, 10 min of inactivation at 95°C, 40 cycles of 15 sec at 95°C, and 45 sec at 60°C. Internal control (EXO), positive, and negative (water) control were included on every run.
Analysis of sequencing data. Reads were processed using a custom bioinformatic pipeline (https://github.com/greninger-lab/covid_swift_pipeline), as described 8 . Briefly, raw FASTQs were trimmed of adapter content and low-quality reads, then aligned to the NCBI reference Wuhan-1 (NC_045512). FASTQs are aligned using BBMap. PCR primers were soft-clipped with Primerclip (https://github.com/ swiftbiosciences/primerclip) and a consensus genome was generated from this alignment using bcftools v1.9 23 . All bam files are manually reviewed to ensure accuracy of alignment parameters, because of the potential for unusually large insertions/deletions that can occur in SARS-CoV-2/coronaviruses. Variants were examined longitudinally with a modified Python script adapted from LAVA (https://github.com/greninger-lab/lava).
A k-mer based approach was used to identify reads corresponding to subgenomic RNAs. After removing adapters and low-quality regions, junctionspanning reads corresponding to individual sgRNAs were filtered using BBDuk (https://sourceforge.net/projects/bbmap/) against a custom FASTA file composed of 30 nucleotides of the 5′ leader sequence and 30 nucleotides of the respective 3' junctions of downstream ORFs 24 , allowing an edit distance of 1. In order to avoid calling of genomic reads, we constrained matches to those that had 10 consecutive matches of 31-mers.
Consensus Orf1a/1b sequences were subject to PAML (version 4.8) positive selection analysis as described 9 .
Steps to produce full-length infectious virus were performed as described previously 11 , with the following minor modifications: Plasmid propagation and digestion were performed as described and excision of DNA bands from the gel was done by using a Blue-Light transilluminator (Accuris). Next, equal molar amounts of the seven plasmids were reduced by half relative to the prior reports to perform the full-length genomic ligation: F1 (0.305 μg), F2 (0.325 μg), F3 (0.375 μg), F4 (0.47 μg) were ligated using T4 DNA ligase with a final volume of 50 μl. Plasmids F5 (0.375 μg), F6 (0.36 μg) and F7 (0.3 μg) were incubated with T4 DNase Ligase in a volume of 40 μl. The primary reactions were incubated at 4°C for 16 h. After this, the two reactions were mixed and volumed up to 100 μl by adding the corresponding additional volume of ligation reagents. The secondary ligation was incubated for 18 h at 4°C. The full-length genomes were purified by using the Genomic DNA Clean and Concentrator kit from Zymo Research according to manufacturer's instructions. Successful ligation products were visualized using a 0.8% agarose gel and standard gel imager.
Next, full-length in vitro RNA transcription was also performed with additional modifications. First, 500 ng of each ligation product was employed. Reactions were performed using the mMESSAGE mMACHINE T7 Transcription Kit (Thermo Scientific) with a cap analog-to-GTP ratio of 1:1. Incubation of the reaction was done at 32°C for 12 h. Following DNAse treatment was performed at 37°C for 15 min. Next, viral RNA was purified by using the RNA Clean and Concentrator kit from Zymo Research according to the manufacturer's instructions. Successful transcription products were visualized using a 0.8% agarose gel and standard gel imager.
RNA transcripts were next electroporated into BHK-N expressing cells (courtesy of Dr. Brett Lindenbach). 15 μg of each mutant transcript was electroporated into BHK cells using the Electro Square Porator Device (BTX). Electroporated cells were seeded into a T75 flask containing 1:2 ratio of BHK-N:Vero-E6 cells. Complete monolayer CPE was evident at 72 h postelectroporation and supernatant was centrifuged and 0.22 um filtered. Vero-E6 were infected with 500 uL of each of the clarified supernatants. Supernatants were collected for each mutant virus at 48 h post infection and tittered as described previously. Infectious clones are available with MTA and appropriate biosafety approvals.
Sequencing and variant calling. RNA was extracted from SARS-CoV-2 mutants using the MagMAX viral/pathogen nucleic acid isolation kit and tested with the multiplexed RT-qPCR variant assay 25 . Using a standard dilution series, Ct values were converted to RNA copies per µL. Extracted RNA was diluted to 1 million copies per µL and used as input into the COVIDSeq Test RUO version for library preparation, in duplicate. Pooled libraries were sequenced on the NovaSeq (pairedend 150) at the Yale Center for Genome Analysis. Single nucleotide variants were called at a minimum frequency threshold of 0.02 and minimum read depth of 400X using iVar (version 1.3.1) 26 , and filtered between duplicate replicates. Introduction of the E802D and E802A was confirmed at a minimum frequency of 0.82-0.9 in initial stocks.
Replication kinetics of icSARS-CoV-2 mutants. Isolates from initial stocks were expanded and used to infect Vero-E6 cells with 0.01 MOI of each E802D and E802A mutant virus to assess their replication kinetics. Supernatants were collected at 8−, 24−, and 48-hours post-infection. Titers for each time point were assessed by conventional plaque assay as described above. Each virus was assessed in triplicate.
Remdesivir resistance assay. Each icSARS-CoV-2 mutant virus (E802D and E802A) was cultured in the presence of increasing concentrations of RDV, ranging from 0.39 to 40 μM. Infections were performed using Vero-E6 in 384-well plates containing the corresponding RDV dilutions in phenol-red free DMEM media with 5% FBS and incubated at 37°C for 4 h. Next, cells were infected with 0.01 MOI of each mutant virus and infected cell frequencies were measured at 48 h postinfection by mNeonGreen expression by high content imaging (Cytation 5, BioTek) 27 . Cell culture supernatants were also collected for plaque assay. Cell viability in uninfected cells was assessed at 72 h post-infection using the CellTiter-Glo kit (Promega) according to the manufacturer's instructions.
Phylogenetic analysis. Phylogenic trees were generated using the Nextstrain software package (version 3.0.6) 28 . 326 background sequences from CT, USA collected April 1st-September 30th 2020 were downloaded from GISAID [www.gisaid.org] (Supplementary Table 3). The tree was rooted using the Wuhan reference strain. Genome sequences which contain the E802D mutation were identified from the GISAID database (accessed on 10/31/21; Supplementary Table 4) and were compared to sequences in the global Auspice dataset.
Flow cytometry. Flow cytometry was performed as described 21  In brief, freshly isolated PBMCs were plated at 1-2 × 10 6 cells per well in a 96-well Ubottom plate. Cells were resuspended in Live/Dead Fixable Aqua (ThermoFisher) for 20 min at 4°C. Following a wash, cells were blocked with Human TruStan FcX (BioLegend) for 10 min at RT. Cocktails of desired staining antibodies were added directly to this mixture for 30 min at RT. For reinfection stains, cells were first washed and supernatant aspirated; then to each cell pellet a cocktail of secondary markers was added for 30 min at 4°C. Prior to analysis, cells were washed and resuspended in 100 μl 4% PFA for 30 min at 4°C. Following this incubation, cells were washed and prepared for analysis on an Attune NXT (ThermoFisher). Data were analysed using FlowJo software version 10.6 software (Tree Star). The specific sets of markers used to identify each subset of cells are summarized in Supplementary Figs. 6 and 8.